#pragma once

#include "OrionFlames.h"
#include "MakeVariation.h"

namespace OrionFlames
{
static const float EPS = (float)1.0e-6;
//EXPORT_STL_VECTOR(OFAPI, Variation*)

MAKE_VARIATION_INVERSE(linear,
	p=t;
,/*inverse*/
	return f(p);
,/* jacobian */
	return 1.0f;
);

MAKE_VARIATION_INVERSE(sinusoidal,
	p=vec2(sin(t.x),sin(t.y));
,/*inverse*/
	if (abs(p.x)>1.0 || abs(p.y)>1.0) return vec4(0.0); /* out of bounds */
	vec4 sum=vec4(0.0);
	vec2 t=vec2(asin(p.x),asin(p.y));
	float img=2; // <- 0-2 is reasonable.  theoretically should be infinity!
	for (float dy=-img;dy<=img;dy++)
	for (float dx=-img;dx<=img;dx++) 
	{
		sum+=f(vec2( t.x+M_PI*2.0*dx,t.y+M_PI*2.0*dy));
		sum+=f(vec2(-t.x+M_PI*(1.0+2.0*dx),t.y+M_PI*2.0*dy));
		sum+=f(vec2(-t.x+M_PI*(1.0+2.0*dx),-t.y+M_PI*(1.0+2.0*dy)));
		sum+=f(vec2( t.x+M_PI*2.0*dx,-t.y+M_PI*(1.0+2.0*dy)));
	}
	return sum;
,/*jacobian*/
	return cos(t.x)*cos(t.y);
);

MAKE_VARIATION_INVERSE(spherical,
	float r2=dot(t,t) + EPS;
	p=t/r2;
,/*inverse*/
	float r2=dot(p,p);
	return f(p/r2);
,/*jacobian*/
	float r2=dot(t,t);
	return 1.0f/(r2*r2);
);

MAKE_VARIATION_INVERSE(swirl,
	float r2=dot(t,t);
	float s=sin(r2); float c=cos(r2);
	p=vec2(
		s*t.x-c*t.y,
		c*t.x+s*t.y
	);
,/*inverse*/
	float r2=dot(p,p);
	if (r2>10.0) return vec4(0.0); /* give up--mipmap can't handle shear... */
	float s=sin(r2); float c=cos(r2);
	return f(vec2(
		s*p.x+c*p.y,
		-c*p.x+s*p.y
	));
,/*jacobian*/
	return 1.0;
);

MAKE_VARIATION_INVERSE(horseshoe,
	float r=(float)(1.0/(length(t)+EPS));
	p=vec2(
		(t.x-t.y)*(t.x+t.y)*r,
		2.0*t.x*t.y*r
	);
,/*inverse: slight numerical problems for big x and small y*/
	float r2 = p.x*p.x + p.y*p.y;
	float r=(float)sqrt(r2);
	float aspect = (p.x + r)/p.y;
	float y = (float)sqrt(0.5*(r2 - p.x*r));
	return f(vec2(y*aspect,y))
	      +f(vec2(-y*aspect,-y));
,/*jacobian*/
	return 2.0f;
);

MAKE_VARIATION_INVERSE(polar,
	p=vec2(
		atan2(t.x,t.y)*M_1_PI,
		length(t)-1.0
	);
,/*inverse*/
	if (p.x>=-1.0 && p.x<=1.0 && p.y>-1.0)
		return f((p.y+1)*vec2(sin(p.x*M_PI),cos(p.x*M_PI)));
	else
		return vec4(0.0);
,/*jacobian*/
	return (float)M_1_PI/length(t);
);

//MAKE_VARIATION(handkerchief,
//	float a=atan2(t.x,t.y);
//	float r=length(t);
//	p=r*vec2(
//		sin(a+r),
//		cos(a-r)
//	);
//,/* jacobian */
//	float r=length(t);
//	return (float)(cos(2.0*r) + 2.0*t.x*t.y/r - r*sin(2.0*r));
//);

MAKE_VARIATION_INVERSE(heart,
	float a=(float)atan2(t.x,t.y);
	float r=(float)length(t);
	//printf("H(%.2f,%.2f -> %.2f,%.2f) ",t.x,t.y,a,r);
	p=vec2(
		r*sin(r*a),
		(-r)*cos(r*a)
	);
,/* inverse */
	float r=(float)length(p);
	float at2=(float)atan2(p.x,-p.y);
	
	/* Figure out the correct multiple(s) of 2PI to add
	 to our arctangent above: */
	float pilo=(float)ceil((-M_PI*r-at2)/(2.0*M_PI));
	float pihi=(float)floor((+M_PI*r-at2)/(2.0*M_PI));
	vec4 sum=vec4(0.0);
	for (float pic=pilo;pic<=pihi;pic++) {
		float a=(float)((at2+pic*2.0*M_PI)/r);
		//printf("HI(%.1f,%.2f,%.2f) ",pic,a,r);
		sum+=f(r*vec2(sin(a),cos(a)));
	}
	return sum;
,/* jacobian */
	return (float)length(t);
);

//MAKE_VARIATION(disc,
//	float a=(float)(atan2(t.x,t.y)*M_1_PI);
//	float r=(float)(length(t)*M_PI);
//	p=a*vec2(
//		sin(r),
//		cos(r)
//	);
///*(self) inverse: This is *only* valid for input points in the 1st quadrant with r<1.0! 
//	float a=length(p)*M_PI;
//	float r=atan2(p.x,p.y)*M_1_PI;
//	return f(vec2(r*sin(a),r*cos(a)));*/
//,/* jacobian */
//	float a=(float)(atan2(t.x,t.y));
//	return (float)(M_1_PI*a/length(t));
//);

//MAKE_VARIATION(spiral,
//	float a=(float)atan2(t.x,t.y);
//	float r=(float)length(t)+EPS;
//	float r1=(float)(1.0/r);
//	p=r1*vec2(
//		cos(a)+sin(r),
//		sin(a)-cos(r)
//	);
//,/* jacobian */
//	float r2=dot(t,t);
//	float r=(float)sqrt(r2);
//	float a=(float)atan2(t.x,t.y);
//	return (float)(1 - r*cos(r - a) + sin(r - a))/r2;
//);

MAKE_VARIATION_INVERSE(hyperbolic,
	/* Was:
	float a=atan2(t.x,t.y);
	float r=length(t) + 1.0e-6;
	p=vec2(sin(a)/r,cos(a)*r);
	*/
	float r2=(float)(dot(t,t)+1.0e-6);
	p=vec2(t.x/r2,t.y); /* simplified w/trig identities */
,/*inverse*/
	float v01 = (float)(1.0/p.x);
	float det = (float)(1.0 - 4.0*p.x*p.x*p.y*p.y);
	if (det>=0.0) {
		float v03 = (float)sqrt(det);
		return f(vec2(((v01*(1.0 - v03))/2.0),p.y))
		      +f(vec2(((v01*(1.0 + v03))/2.0),p.y));
	} else { /* negative determinant--no solutions */
		return vec4(0.0);
	}

,/* jacobian */
	/* Was:
		float a=atan2(x,y);
		return cos(2*a)/r2;
	
	I'm sure there's some way to simplify this, e.g.,
	with cos(2*a) = 2 (cos a)^2-1
	cos(a)*r=y; cos(a)=y/r; 
	cos(2a)=2 y^2/r^2-1
	*/
	float r2=dot(t,t);
	return (float)((1.0-2.0*t.y*t.y/r2)/r2);
);

//MAKE_VARIATION(diamond,  /* has an inverse, but like 32 roots! */
//	float a=atan2(t.x,t.y);
//	float r=length(t);
//	p=vec2(sin(a)*cos(r),cos(a)*sin(r));
//,/* jacobian */
//	float r2=dot(t,t);
//	float r=sqrt(r2);
//	return (float)((cos(2.0*r)+2*t.y*t.y/r2-1)/(2.0*r));
//);

//MAKE_VARIATION(ex,
//	float a=(float)atan2(t.x,t.y);
//	float r=(float)length(t);
//	float n0=(float)sin(a+r);
//	float n1=(float)cos(a-r);
//	float m0=n0*n0*n0*r;
//	float m1=n1*n1*n1*r;
//	p=vec2(m0+m1,m0-m1);
//,/* jacobian */
//	float r2=dot(t,t);
//	float r=(float)sqrt(r2);
//	float s2r=(float)sin(2.0*r);
//	float c2r=(float)cos(2.0*r);
//	float inner=s2r+t.x*t.y/r2;
//	return (float)(1.5/r*(6.0*t.x*t.y+r*c2r-3*r2*s2r)*inner*inner);
//);

/* FIXME: julia is random, but only two possible outcomes--
  so split incoming maps, make other map have negative wval 
  to flip output around */
MAKE_VARIATION_INVERSE(julia,
	/* Trig version:
	float a=atan2(t.x,t.y);
	float r=length(t);
	p=sqrt(r)*vec2(cos(a*0.5),sin(a*0.5));
	*/
	
	// Man, trig identities are friggin' weird.  Half-angle gives:
	float r=length(t);
	p=vec2(sqrt(0.5*(r+t.y)),sqrt(0.5*(r-t.y)));
	if (t.x<0.0) p.y=-p.y; /* flip sign */
	
,/* inverse */
	if (p.x<0) return vec4(0.0); /* that point's not one of ours */
	/* Triggy inverse is this:
	float a=atan2(p.y,p.x);
	float r2=dot(p,p);
	return f(r2*vec2(sin(a*2.0),cos(a*2.0)));
	*/
	// Mathematica claims this is also the inverse:
	return f(vec2( 2*p.x*p.y, (p.x - p.y)*(p.x + p.y) ));
	
,/* jacobian */
	return (1.0f/(4.0f*length(t)));
);

MAKE_VARIATION_INVERSE(bent,
	vec2 n=t;
	if (n.x<0.0) n.x*=2.0;
	if (n.y<0.0) n.y*=0.5;
	p=n;
,/* inverse */
	vec2 n=p;
	if (n.x<0.0) n.x*=0.5;
	if (n.y<0.0) n.y*=2.0;
	return f(n);
,/* jacobian */
	float scale=1.0;
	if (t.x<0.0) scale=2.0;
	if (t.y<0.0) scale*=0.5;
	return scale;
);

/* FIXME: waves somehow adjusts parameters based on coeffs */
MAKE_VARIATION_INVERSE(fisheye,
	float r=(float)length(t);
	r=2/(r+1);
	p=r*vec2(t.y,t.x); //<- note backwards!
	/* length (p) = s
		s=r*2/(r+1)  (w/ extra copy of length(t) from vec2)
		sr+s=r*2   (solve for r)
		r*(s-2)=-s
		r=s/(2-s)
	*/
,/*inverse*/
	float s=(float)length(p);
	float r=(float)(1.0/(2.0-s));
	if (r>=0.0) /* plausible fisheye output: */
		return f(r*vec2(p.y,p.x));
	else /* outside fisheye region */
		return vec4(0.0);
,/*jacobian*/
	float r=(float)length(t);
	float r1=(float)(1.0+r);
	return (float)(4.0/(r1*r1*r1));
);

/* FIXME: popcorn needs coeffs too */
MAKE_VARIATION_INVERSE(exponential, 
	float dx=(float)exp(t.x-1.0);
	float dy=(float)(M_PI*t.y);
	p=dx*vec2(cos(dy),sin(dy));
,/*inverse*/
	float r=(float)length(p);
	float a=(float)atan2(p.y,p.x);
	return f(vec2(log(r)+1.0,a*M_1_PI));
,/*jacobian*/
	return (float)(M_PI*exp(2.0*t.x-2.0));
);

MAKE_VARIATION_INVERSE(power,
	float a=(float)atan2(t.x,t.y);
	p=vec2(cos(a),sin(a))*pow(length(t),sin(a));
,/*inverse*/
	float a=(float)atan2(p.y,p.x);
	float r=pow((float)length(p),(float)1.0/sin(a));
	if (r==0.0 || r+1.0==r) /* r has exploded--give up */
		return vec4(0.0);
	else /* returned a reasonable value */
		return f(r*vec2(sin(a),cos(a)));
,/*jacobian*/
	float ir=(float)(1.0/length(t));
	float r=(float)(1.0/ir);
	return (float)(pow((float)r,(float)(2.0*(t.x*ir)-2.0))*(t.x*ir));
);

//MAKE_VARIATION(cosine,  /* has 16 families of solutions */
//	float n=(float)(t.x*M_PI);
//	p=vec2(cos(n)*cosh(t.y),-sin(n)*sinh(t.y));
//,/*jacobian*/
//	return (float)(0.5*M_PI*(-cos(2.0*M_PI*t.x)+cosh(2.0*t.y)));
//);
/* FIXME: rings needs coeffs */

class OFAPI VariationList
{
public:
	VariationList()
	{
		ClearVariations();
		m_VariationVector.push_back(new Variation_linear());
		m_VariationVector.push_back(new Variation_sinusoidal());
		m_VariationVector.push_back(new Variation_spherical());
		m_VariationVector.push_back(new Variation_swirl());
		m_VariationVector.push_back(new Variation_horseshoe());
		m_VariationVector.push_back(new Variation_polar());
		//m_VariationVector.push_back(new Variation_handkerchief());//forward only?
		m_VariationVector.push_back(new Variation_heart());
		//m_VariationVector.push_back(new Variation_disc());//forward only?
		//m_VariationVector.push_back(new Variation_spiral());//forward only?
		m_VariationVector.push_back(new Variation_hyperbolic());
		//m_VariationVector.push_back(new Variation_diamond());//forward only?
		//m_VariationVector.push_back(new Variation_ex());//forward only?
		m_VariationVector.push_back(new Variation_julia());
		m_VariationVector.push_back(new Variation_bent());
		m_VariationVector.push_back(new Variation_fisheye());
		m_VariationVector.push_back(new Variation_exponential());
		m_VariationVector.push_back(new Variation_power());
		//m_VariationVector.push_back(new Variation_cosine());//forward only?
		m_VariationVector.push_back(NULL);
	}

	~VariationList()
	{
		ClearVariations();
	}

	static unsigned int Size() { return m_VariationVector.size(); }
	static Variation* Get(int index) { return m_VariationVector[index]; }

	static Variation* LookupVariation(string name)
	{
		for (unsigned int i = 0; i < m_VariationVector.size() && m_VariationVector[i] != NULL; i++)
			if (name == m_VariationVector[i]->getName())
				return m_VariationVector[i];
		
		throw ("LookupVariation() : Bad Variation lookup " + name);
		return NULL;
	}

	static int LookupVariationIndex(string name)
	{
		for (unsigned int i = 0; i < m_VariationVector.size() && m_VariationVector[i] != NULL; i++)
			if (name == m_VariationVector[i]->getName())
				return i;
		
		throw ("LookupVariationIndex() : Bad Variation lookup " + name);
		return -1;
	}

	static int LookupVariationIndex(Variation* var)
	{
		for (unsigned int i = 0; i < m_VariationVector.size() && m_VariationVector[i] != NULL; i++)
			if (var == m_VariationVector[i])
				return i;
		
		throw ("LookupVariationIndex() : Bad Variation lookup " + var->getName());
		return -1;
	}

private:
	static void ClearVariations()
	{
		for(unsigned int i = 0; i < m_VariationVector.size(); i++)
		{
			if (m_VariationVector[i] != NULL)
				delete m_VariationVector[i];

			m_VariationVector[i] = NULL;
		}

		m_VariationVector.clear();
	}

	static vector<Variation*> m_VariationVector;
};
}